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We present a Monte Carlo simulation of an effective theory for local 
Polyakov loops at finite temperature and density. The sign problem is 
overcome by mapping the partition sum to a flux representation. We de- 
termine the phase diagram of the model as a function of the temperature 
and the chemical potential. 
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1. Introduction 

Lattice QCD is a powerful tool to address non-perturbative phenomena 
quantitatively and in principle is one of the most appropriate techniques to 
explore the QCD phase diagram. However, at finite chemical potential the 
fermion determinant becomes complex and it can not be used as a Boltz- 
mann weight in Monte Carlo simulations. Alternative approaches such as 
reweighting, power series expansion, strong coupling/large mass expansion 
or analytic continuation from imaginary chemical work only for small chem- 
ical potential leaving the rest of the phase diagram unexplored. For true 
progress with QCD thermodynamics on the lattice new ideas are necessary. 

In this article we explore the phase diagram of the SU(3) spin model pQ, 
where the degrees of freedom are traced SU(3) valued spins (local Polyakov 
loops) as a function of temperature and chemical potential. This effective 
theory can be derived from full QCD using strong coupling expansion for 
the gluon action and hopping expansion for the fermion determinant. It 
is motivated by the relation of the deconfinement transition and center 
symmetry of pure gauge theory [2] . From the fermion determinant one takes 
into account a center symmetry breaking term which couples the chemical 
potential /i and gives rise to a sign problem at finite /i. However, in this 
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model the sign problem can be solved by exactly rewriting the partition 
sum in terms of flux variables [H [5] . 

2. Center effective theory 

The action of the center effective theory has the form 




P(x)P{x + z>)* + c.c. + k e^P{x) + e-»P(x)* . (1) 



The degrees of freedom P(x) are the traced SU(3) variables P(x) = Tr 
L{x) with L(a?) € SU(3) attached to the sites x of a three-dimensional cubic 
lattice with periodic boundary conditions. By v we denote the unit vector 
in ^-direction, with z/ = 1, 2, 3. The first term of the action, i.e., the nearest 
neighbor interaction term, can be obtained as the leading contribution in 
the strong coupling expansion of the gauge action. This term is invariant 
under center transformations P(x) — > zP(x) with z £ Z3. The parameter r 
depends on the temperature (it increases with T) and is real and positive. 
The second term, referred to as the magnetic term, is obtained as the leading 
contribution in the hopping expansion (large mass expansion) of the fermion 
determinant. The real and positive parameter k is proportional to the 
number of flavors and depends on the fermion mass (it decreases with m q ). 
The magnetic term breaks center symmetry explicitly and is complex when 
the chemical potential \i is non-zero, thus generating a sign problem. 

The grand canonical partition function of the model described by ([T]) is 
obtained by integrating the Boltzmann factor e~ s ^ over all configurations 
of the Polyakov loop variables. The corresponding measure is a product 
over the reduced Haar measures dL(x) at the sites x. Thus 

Z = TT f dL(x)e~ s ^ = [ D[L]e- s ^ . (2) 
j. JSU(3) J 

Equations ([T]) and t^fy define the SU(3) effective theory. 



3. Solving the sign problem 

To overcome the sign problem we apply high temperature expansion tech- 
niques and map the theory onto a flux representation, where the partition 
function is rewritten in terms of new degrees of freedom, so called flux vari- 
ables. Here we outline the general strategy for the derivation of the flux 
representation (for the details see [5]). The general steps are: 
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1. Write the Boltzmann weight in a factorized form and expand the ex- 
ponentials for individual links and sites. 

• Nearest neighbor term (links): 

e rP( X rP( X+ u) ^ Y^i^[P{xfP{x + i>)\ 1 "" 

- lx t V' 

• Magnetic term (sites), we use r] = ne^ 1 and rj = ne~^: 

Sx- — S x - 

Sx Sx 

2. Rewrite the partition function as: 

_where f(x) = YjI=i[ 1 x,v + l x-u,u] + s x and f(x) = Yl=i\ l x-o,u + 
lx,u] + s x denote two types of fluxes at a site x of the lattice. 

3. After integrating out the SU(3) variables L{x) |3J, the new form of 
the partition sum depends only on the flux variables: 

• Dimers l x ,uJx,u £ [0, +oo[ , living on the links (x, v). 

• Monomers s x ,s x € [0, +oo[ , living on the sites x. 

4. The flux variables l x ,uilx,ui s x ,s x are the new degrees of freedom and 
Ylun S{ s s} denotes the sum over all their configurations. The flux 
variables are subject to a constraint which forces the total flux f(x) — 
f(x) to be a multiple of 3 at each site x. 

4. Numerical analysis 

For the analysis we performed simulations with a local Monte Carlo update 
on 10 3 , 16 3 and 20 3 lattices and focused on the bulk observables internal en- 
ergy U and the magnetization P (which is identified with the Polyakov loop 
of QCD), as well as their fluctuations C (heat capacity) and xp (Polyakov 
loop susceptibility). 

First we performed several checks of the flux representation and the algo- 
rithm. Fig. [T] shows that for small [i and r the data obtained from the 
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Fig. 1. P (left) and X p (right) for n = 0.1,0.04 and 0.005 and r = 0.001. We 
compare the results from the Monte Carlo simulation on a 10 3 lattice (circles) and 
perturbative expansion in r (lines). 



simulation (circles) nicely approaches the analytical results from a pertur- 
bative expansion in r (lines). We plot P and xp f° r r = 0.001 and three 
different values of function of the chemical potential. The same com- 



parison is shown in Fig. 3(a) , where the solid curves at the bottom are the 



positions of the maxima from the perturbative expansion. 

We also compared our results to other approaches. Fig.[2]shows that the flux 
results and the data from a complex Langevin calculation |6j agree very well, 
and for fi = also with the results from a conventional simulation in the spin 
representation. The discrepancy at r = 0.132 is solved when a higher-order 
algorithm is used for the two values, fi 2 = 0.0 and 0.2 (crosses). 

To explore the phase boundaries in the r-/x plane, we identified the posi- 
tions of the maxima of xp an d C. Subsequently we used two methods to 
determine the nature of the transitions: first we studied the histograms of 
U and P to search for a double peak behavior characteristic of a first or- 
der transition, and secondly we analyzed the volume scaling of C and xp- 



Fig. 3(a) shows the positions of the maxima of xp m the t-\i plane for 
k = 0.1, 0.04, 0.02 and 0.005. We find that there is a first order phase 
transition for small fi and K < k c (triangles), while the rest of the transi- 
tion lines shows a crossover behavior (circles). Our estimate for the critical 
point for /z = is (t c ,k c ) = (0.1331(1), 0.016(2)). This value is different 
from a mean field analysis of the SU(3) spin model [7\ where the critical 
point was reported to be at k = 0.059. However, in [8] it was shown that 
when considering higher order corrections the value k c from the mean field 
approach decreases. Fig. |3(b)| shows the positions of the maxima of xp 
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Fig. 2. Comparison of (P + P*) /2V from the flux simulation (filled symbols) to the 
results from the complex Langevin approach (empty symbols and two high accuracy 
data points are marked with crosses). For \l = we also added the results from 
a simulation in the conventional spin approach (asterisks). We compare data at 
different values of r as a function of /j, 2 for k = 0.02 on lattices of size 10 3 . 
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Fig. 3. (a) Left: Phase diagram obtained from the maxima of xp f° r 4 values of k. 
The horizontal line marks the critical t for k = 0, and the curves at the bottom are 
the results from a r expansion. The red point is the critical end point for k = 0. 
(b) Right: Comparison of the phase boundaries obtained from the maxima of xp 
and C for three values of k. 



and C, demonstrating that the crossover region (manifest also in different 
positions for the maxima of xp anu C) becomes wider with increasing [i. 
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5. Conclusions 

We have studied an effective theory for the Polyakov loop at finite temper- 
ature and density. Mapping the theory onto a flux representation enables 
us to have a model free of the sign problem and opens the possibility to use 
Monte Carlo techniques. For large values of n (physical case) the transition 
is of a smooth crossover type and we conclude that center symmetry alone 
does not provide a mechanism for first order behavior in the QCD phase 
diagram. 

We also compared the recently published results from a complex Langevin 
simulation of the SU(3) spin model [6] to the data from our flux simulation. 
We find very good agreement between the two methods which is a valuable 
test for both, the flux and the complex Langevin approach. 
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